% Monte Carlo simulation:
%
% Joint hypothesis test: alpha_1 & beta_1 of GECM are significant
%               or
% LRM is significant
% 
% Estimating GECM with either fractionally integrated data: (d(y) = d(x))
%               or 
% GECM with near-integrated data
% 
% Requires fracdiff.m - available from Katsumi Shimotsu
% http://shimotsu.web.fc2.com/codes/elwcode10.zip (zip file)
% 
% Requires Ols.m - found in the Oxford MFE Toolbox, available
% from Kevin Sheppard
% http://www.kevinsheppard.com/MFE_Toolbox 
%____________________________________________________________

clear;

tic 

% Set Seed % 
%randn('state',1234);       % for exact replication of NI series
randn('state',538151819);   % for exact replication of FI series

ns = [100 200 300 400 500];  % number of observations
d1s = [.2];  % d -- after fractional differencing, series is at this level of d 
rep = 1000; 
phix = [.9];
phiy = [.9];

nd = length(d1s); 

for i = 1:nd
    d1 = d1s(i);
    
    dhat = [];
    summary = [];
    
    for j=1:length(ns);
        n = ns(j);
        
        for k = 1:rep
            
            u1 = randn(n+50,1); % random draw with extra 50 observations
            u2 = randn(n+50,1);
            
            % Simulate Process - Choose one (FI or NI), comment out other %
            
            % Near Integrated (NI) - add AR parameter to the series (1,0,0)
%             ar_x = filter(1,[1 -phix],u1);
%             ar_y = filter(1,[1,-phiy],u2);
%             % Copy each series %
%             x = ar_x;
%             y = ar_y;
            
            % Fractionally Integrated (FI) - (0,d,0)
            fi_x = fracdiff(u1,-d1);
            fi_y = fracdiff(u2,-d1);
            % Copy each series %
            x = fi_x;
            y = fi_y;
            
            % drop the first 50 observations
            x(1:50,:) = [];    
            y(1:50,:) = [];
            
            % Generate lags for GECM
            lagY = [y(1:end-1)];
            lagX = [x(1:end-1)];
            
            % Create matrix for x and difference DV
            x1 = [lagY, diff(x), lagX]; 
            d_y = diff(y);
            
            % Naive Regression (no accounting for FI)
            [b,~,~,VCV] = ols(d_y,x1,1);   % GECM w/o accounting for FI
            tstat = (b./sqrt(diag(VCV))); % tstat for plain SEs
            % tstat order = [constant, lagY, diff(x), lagX]
            
            % Estimate LRM, SE of LRM and t-stat of LRM (using formula from D&K)
            LRM = b(4)/b(2); % lrm is (lagX / lagY)
            LRMse = sqrt((1/b(2)^2)*VCV(4,4) + (b(4)^2/b(2)^4)*VCV(2,2) - 2*(b(4)/b(2)^3)*VCV(4,2));
            LRMt1 = LRM/LRMse;
            
            % T-statistics
            tval(k,:) = tstat(2);   % tstat - lagY
            tval2(k,:) = tstat(4);  % tstat - lagX
            tval_lrm(k,:) = LRMt1;  % tstat - LRM
            
            % Significance
            tsig = (abs(tval) > 1.645 & abs(tval2) > 1.96)+0;    % +0 to convert logical to double
            tsig_lrm = (abs(tval_lrm) > 1.645)+0;
            
            sigi1 = tabulate(tsig);
            sig_lrm = tabulate(tsig_lrm);

        end
        
        % % % % % % % % % % % % % % % % % % % % % % %
        
        % Compute percentage of significant ECM and LRM

        sigpct_ECM = 1 - (sigi1(1,2)/rep);
        sigpct_LRM = 1 - (sig_lrm(1,2)/rep);

        % Summarize and display the results (Comment out one not being used)
        
        % FI %
         fprintf(1, '\t d \t %%ECM \t %%LRMsig   N  \n ')
         summary = [summary; d1' sigpct_ECM' sigpct_LRM' n'];
        
        % NI % 
%            fprintf(1, '\t phix \t phiy \t%%ECMsig  %%LRMsig   N  \n ')
%            summary = [summary; phix' phiy' sigpct_ECM' sigpct_LRM' n'];
          
         disp(summary); 
        
    end

end
toc
